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Abstract 

In simulating continuum model fluids that undergo phase separation and criticality, 
significant gains in computational efficiency may be had by confining the particles to 
the sites of a lattice of sufficiently fine spacing, ao (relative to the particle size, say 
a). But a cardinal question, investigated here, then arises, namely: How does the 
choice of the lattice discretization parameter, (" = a/oo, affect the values of interesting 
parameters, specifically, critical temperature and density, Tc and pc? Indeed, for small 
C 4-8) the underlying lattice can strongly infiuence the thermodynamic properties. 
A heuristic argument, essentially exact in d = 1 and d = 2 dimensions, indicates 
that for models with hard-core potentials, both Tc(C) and Pc(C) should converge to 
their continuum limits as for d < 3 when C oo; but the behavior of the 

error is highly erratic for d > 2. For smoother interaction potentials, the convergence 
is faster. Exact results for d = 1 models of van der Waals character confirm this; 
however, an optimal choice of C can improve the rate of convergence by a factor 
For d > 2 models, the convergence of the second virial coefficients to their continuum 
limits likewise exhibit erratic behavior which is seen to transfer similarly to Tc and pc', 
but this can be used in various ways to enhance convergence and improve extrapolation 
to C = oo as is illustrated using data for the restricted primitive model electrolyte. 



I Introduction and General Considerations 



Understanding the thermodynamic behavior of complex fluids, such as electrolytes, poly- 
mers, colloids, etc., has been a key issue in soft-condensed-matter physics in recent years. 
Simulating realistic continuum models of such fluids with the aid of powerful computers can 
be informative; however, very large computing resources may be required. Even determining 
accurately the critical parameters of the simplest model of an electrolyte, namely the re- 
stricted primitive model (RPM), has been an outstanding problem for more than a decade.^ 
To ease the difficulties of simulation, various computer algorithms and special techniques 
have been developed. Particularly helpful is the fine-lattice discretization approach intro- 
duced by Panagiotopoulos and Kumar^ that has provided an efficient way of simulating the 
thermodynamic properties near phase separation and criticality in various model fluids. 
In this method, as illustrated in Figure 1, the center of any particle may reside only on the 
sites of a lattice of sufficiently fine spacing qq, relative to the particle size, a, which it is 
usually convenient to take as the diameter of the intermolecular repulsive core. A crucial 
parameter is then the lattice-discretization ratio which we define generally by 

C = a/ao > 1, (1) 

where, as is worth stressing at the outset, C need not be an integer. Systems with C = 1 
correspond to the simplest lattice gas models with, typically, single-site hard cores that 
merely prevent double occupancy of a lattice site. On the other hand, when ^ — > oo one 
clearly recaptures the full continuum version of the model fiuid under study. In practice one 
observes that even for moderately low values of (, say in the range 5-10, the properties of 
the discretized systems are close to those of their continuum limits.^'^ 

The advantage of simulating such discretized systems rather than their continuum coun- 
terparts is that one can readily generate a look-up table of the interactions between two 
particles at all allowable distances (up to a maximum of order the system size, L) prior 
to undertaking the simulations. This then speeds up the computations markedly: thus for 
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the discretized RPM one gains a factor of speed as large as 100 relative to the continuum 
model. ^ On the other hand, one drawback of the approach is that one cannot simulate using 
arbitrarily large values of ( owing to the limitations on the capacity of fast-access memory 
in which to store the interaction data. Hence, some compromise is inescapable — especially 
when large system-sizes, L, are demanded^" — in order to select a reasonably low value of ( 
that will still ensure satisfactory estimates of continuum properties. 

To this end we hence ask: How does the lattice-discretization parameter, affect the 
thermodynamics of the model fluids under study? Of particular interest, and a focus of this 
article, is the ^-dependence of the critical parameters of systems that undergo phase sepa- 
ration and gas-liquid or fluid-fluid criticality. In particular, how do the critical temperature 
and density, and Pc, depend on (1 More specifically, one may expect that Tc{C) and 
Pc(C) will converge to their continuum limits as and 1/C^'', respectively, with conver- 

gence exponents ipT and ipp that arc clearly of interest. It is, indeed, reasonable to expect 
i^T — i^p — for general cases; but, how does -0 depend on the details of the models? Fur- 
thermore, understanding these issues may provide a route to reliably estimating properties 
such as Tc and precisely via an extrapolation to the continuum limit {C, oo) based on 
simulations at moderately low values of C,. 

Numerical studies of the (^-dependence of Tc and pc have previously been performed for 
certain model fluids. "^'^'^ For the RPM — hard-core ions carrying charges ±go — Panagiotopoulos'' 
concluded that both the critical temperature and density approach their continuum limits 
approximately as 1/^^ for values of C, from 5 up to 15. On the other hand, uncharged models 
with strong but smoothly rising, e~"^ repulsive interactions have exhibited faster rates of 
convergence of T'c(C); namely, the deviations T'c(C) — 31 (oo) appeared to scale as l/C^ with 
■0 = 6 ± 2 for C, in the range (5, 10).^ [Moreover, for the models examined the estimates of 
Pc(C) displayed rather small dependence on C,.\ In other calculations,^ it was observed that 
a charged dumbbell model with hard-core repulsions exhibited the same convergence rate 
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as the RPM, namely with an exponent ~ 2. Later, rather more precise simulations^ con- 
firmed the approximate overall convergence of T^. and Pc for the RPM but the detailed 
variation with ( was, in fact, found to be quite erratic or 'noisy' at a level significantly larger 
than the simulation uncertainties. How can one understand these numerical observations 
from an analytical or general theoretical viewpoint? 

To answer this question, it is essential to understand the nature of the truncation errors 
of the thermodynamic and statistical properties of the discretized system relative to its 
continuum counterpart as ( changes. More concretely, in the grand canonical ensemble, 
for example, the partition function, E{T,ii), can be obtained by integrating the Boltzmann 
factor for the A^^-particle Hamiltonian Hat over the full phase space and by summing on 
A^. On the other hand, in the discretized system, the configuration space is limited to the 
lattice points of spacing qq — a/C; thus the spatial integrations defining S must be replaced 
by appropriate summations over lattice configurations, yielding a S''(T, /x) which deviates 
from the continuum limit, H°°(T, //). The grand thermodynamic potential from which the 
thermodynamic quantities can be derived is then proportional to InS. How does C affect the 
partition function (and thus the corresponding thermodynamics) for the discretized system? 

To simphfy the situation and to gain a first insight into the character of the C-convergence 
of Tc and pc, it is instructive, following ref 8, to examine an approximate, van der Waals 
(vdW) description of the equation of state for a model fiuid. Accordingly, for a d-dimensional 
fiuid with pair interaction potential (p{r) — (fio{r) + (fii{r), where (po is repulsive and of short 
range, while ipi is smooth, attractive and long ranged, let us consider the vdW-type equation 
of state. 



which will be valid semiquantitatively and, in general, exhibit (classical) criticality. In this 
equation Bq{T) and Bi{T) are partial second virial coefficients given by 



p/phsT = [1 - pBo{T)]-' + pB,{T), 



(2) 




(3) 
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while for the full second virial coefficient, B'^iT) — Bo{T) + Bi{T), one simply replaces (pi{r) 
by the total potential (p{r). For a discretized system with ( = a/ao, the B^ integrals must 
be replaced by sums, 



over integral lattice vectors, n, and hkewise for B^{T). The critical temperature and density, 
Tc(C) and pdC), are then readily derived from Bq{T) and B^{T). Thus it is straightforward to 
see that the convergence of the B^{T) to their continuum limits, Bf°{T) — characterized, say, 
by an exponent ■0s — rather directly governs the behavior of T'c(C) and Pc(C)- particular, 
in the absence of special or "accidental" cancellations, one will have 4>t = '4^p = 4>b- 

Of course, we do not expect this approach to yield quantitatively accurate information 
since the true critical behavior in two- and three-dimensional systems deviates from the clas- 
sical, vdW paradigm. However, we do expect the results to be correct at a semiquantitative 
level as, indeed, are the estimates of Tc and pc obtained from the van der Waals equation 
itself.^ Similarly, we anticipate, and will demonstrate numerically and analytically, that the 
character of the asymptotic behavior as embodied in the exponents V', etc., is in general 
preserved. 

Nevertheless, it is not unreasonable to question whether the existence of the nonclassical 
critical exponents /?<|,7>l,i/>|, etc., for the true thermodynamic coexistence curve, 
compressibility, and correlation length, etc., might not change the convergence exponents ■0 
for Tc(^) and Pc(C)- This is an interesting question which we address in Appendix A. The 
arguments there indicate that in leading order when C, increases, the behavior of the critical 
parameters will be dominated by the convergence of low-order integrals over Boltzmann 
factors, as epitomized in the virial coefficients, and thus will in be insensitive generically to 
the nonclassical subtleties of the true critical behavior. 

To follow the hne of investigation proposed we must thus ask: How rapidly will the B^{T) 




(4) 



n 
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converge to their limits? Specifically, if we define the relative truncation error^ via 

E.(C;T) = 4(T)/5r(T)-l, (5) 

and consider the approximation of the integrals by the sums it is rather clear that 
the rate of decay of -Ej(C) to as C increases will be determined by the continuity and 
smoothness of the integrands fi{r). In particular, if fi{r) or any of its low-order derivatives 
exhibit discontinuities, we should expect these to play a dominant role. It follows at once 
from our description of foir) and ifi{r) that the strongly varying repulsive potentials and, 
thence, Eo{(), should be the major players. Indeed, some of the explicit calculations we 
present below will amply illustrate this inference. Already we can see why the RPM and 
the hard-core dumbbell models should, by virtue of the corresponding discontinuities in the 
relevant Mayer /(r) functions, have been expected to exhibit a slower convergence than the 
uncharged models with smooth repulsive potentials. 

Because they lead to the slowest convergence and because of their ubiquity in theoretical 
modeling (albeit primarily for analytical and conceptual convenience) let us focus on the 
effects of hard-core repulsive potentials. A few moments thought shows that 2i?Q(T)/ag is 
just the number of lattice points, say NalC), contained within a c?-dimensional sphere centered 
on a lattice site and of radius ( (using units of the lattice spacing ao): see Figure 1. Likewise 
2B^(T)/aQ is the volume of the corresponding sphere, say Vd{C). If, as for simplicity we will 
suppose henceforth, attention is confined to square, simple cubic, and hypercubic lattices 
[with lattice vectors n = (jir) (r = 1,2, ■ ■ ■ ,d) and Ur = 0, ±1, ±2, ■ ■ ■], one soon realizes 
that the behavior of Nd{C) as ( increases is a classic problem in the geometry of numbers. 

A heuristic argument, presented in Appendix B for the insight we feel it conveys, suggests 
that the truncation error -E'o(C) decays as for d <3, while for > 3 its decay follows 

a power law, independent of d. This argument can be checked immediately for = 1 
since for any non-integral ( in a linear lattice we have |Vi(C) — Ni{()\ < 2 while Vi{() = 2( 
so that -E'o(C) decays as On the other hand, for c? > 4 it has been proven rigorously^^ 
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that -E'o(C) decays as whatever the value of d}'^'^^ For d — 2, Bleher, with Dyson and 
Lebowitz/^'-'^^ has cstabhshed that Eq{Q is of magnitude but has an erratic or noisy 

behavior. Finally, Bleher and Dyson^^'^^ have shown that o? = 3 is a borderline dimension at 
which -E'o(C) contains a logarithmic factor and, indeed, varies in magnitude as (lnC)''^^^/C^- 
These results confirm the heuristic arguments — see section II and Appendix B — except for 
the appearance of the factor \/\nX, when d = 3: this presents an interesting analogy with the 
theory of critical phenomena where the free energy and other properties invariably contain 
logarithmic factors at borderline dimensions. 

Indeed, as we have indicated, the convergence of Tc{C,) and Pc(C) to their limits observed 
numerically in the RPM and the hard-core dumbbell model agrees with our analyses. How- 
ever, the logarithmic correction factor is too difficult to detect, in light of the limited (small 
integer) values of C examined. The faster convergence of T'c(C) and Pc(C) iii the uncharged 
models^ may clearly be attributed to the smoothness of the repulsive potentials adopted. 
On general grounds, well illustrated in section III for [d — l)-dimensional systems, one ex- 
pects that as the integrands /o(r) and /i(r) exhibit discontinuities in successively higher 
derivatives 0, 1, 2, • • • , so the truncation errors will decrease by further factors of l/Q. Thus 
for Boltzmann factors ex^[—ip{r)/k-QT] with continuous second derivatives one might ex- 
pect truncation errors decreasing faster than for d = 3 dimensions. This argument 
must be viewed with caution, however, even as regards potentials with many or all deriva- 
tives continuous: indeed, it is most probable that the observed exponent estimates for these 
smooth-potential models, namely ip = 6±2, represent effective (rather than truly asymptotic) 
exponents that are relevant to the finite range of ( explored and to the relative sharpness of 
the Boltzmann-factor variation: see section IV. 

The balance of this article is organized as follows. In section II the behavior of the partial 
second virial coefficient Bq for a hard-core potential is discussed in further detail focusing, 
in particular, on d — 2 and 3. In section III, following Fisher and Widom,^^ exact results for 
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the critical parameters of one-dimensional discretized van der Waals models are obtained for 
various short-range potentials. The convergence of the calculated values of Tc{() and Pc(C) 
demonstrates in explicit fashion the role of discontinuities and the smoothness of the pair 
potentials. The truncation errors in two and three dimensions are considered in section IV 
with illustrations based, in particular, on simulation data for the RPM.® A brief summary 
is presented in section V. 

II Hard-core Second Virial Coefficient 

For subsequent applications to extrapolating to the continuum limit the critical points of 
the discretized RPM, or other model fluids with hard-core interactions, we consider in this 
section the hard-core part of the second virial coefficient, namely, Bq. As pointed out, this 
is equivalent, via suitable normalization, to determining the number of lattice points Nd{C) 
inside a ci- dimensional sphere of radius R = (gq. When C ^ cxd, it is obvious that Nd{C) 
increases like the reduced volume, Vrf(C) ~ C^, of the sphere; but the crucial question is: How 
does Nd{C) approach the continuum behavior? One may write 

iVd(C) = Ki(C) + MC), (6) 

where the deviation nd{C) satisfies nd{C)/C'^ ^ as C — ^ oo. More specifically we wish to 
know how nd{C) depends on (. 

For d=l, it is trivial that nrf~0(l) so that Eo{() =nd{()/Vd{() defined as in eq El 
decays as see Figure 2(a) which displays a scaled plot. Note, first, that -Eo(C) is a 
piecewise continuous function, locally always decreasing but with infinitely many positive 
jump discontinuities. These features, in fact, characterize the variation of -Eo(C) fo^' all d,: see 
Figure 2(b). Second, observe that Eq{() vanishes when ( takes half-integer values: this, of 
course, is special to d=l. In the next section these values for ( will be shown to provide faster 
convergence for the critical parameters of discretized one- dimensional models. A selection 
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of similar special ( values for d — 2 and d — 3 dimensions is presented in Table I below. 

To understand the convergence of the truncation error Eo{() for d > 1, consider a d- 
dimensional ball of radius ( centered at the origin site of a simple hypercubic lattice. The 
error term, ^^(C), then arises near the surface of the ball which cuts elementary lattice cells 
into two parts. A heuristic argument presented in Appendix B supposes that the cutting 
of nearby cells is only weakly correlated: in that case the root mean square magnitude of 
UdiC) grows like (('^-^y^. This imphes that Eo{() should decay as This argument, 

however, is inadequate for o? > 3 where the exact results show that £^o(C) decays only as 
with a convergence exponent independent of d as recounted in section I.^^ In other 
words, the error function UdiC) grows as C^"^ for d > 3, i.e., much faster than C^*^"^)/^ as 
suggested by the heuristic argument with limited correlations between nearby cells. Indeed, 
as demonstrated in the second part of Appendix B, the reason for the distinct behavior is 
that the surface of a sphere is effectively flatter in higher dimensions {d > 3) so that the 
lattice cells cut by the sphere are strongly correlated over regions of the sphere having d — 2 
dimensions of order the radius R. 

For our present purposes the most interesting cases are d = 2 and 3, as appropriate for 
reahstic model fluids. For d — 2, Bleher, Dyson and Lebowitz^^ demonstrated some time 
ago, in a completely different physical context, that n2{() fluctuates rather wildly growing in 
magnitude as C^/^. The truncation error, £^o(C)) then decays as as illustrated in Figure 

2(b). This is, indeed, in accord with our heuristic argument. However, Bleher et al}^ went 
further by showing that the scaled fluctuation, 'n-2(C)/C^^^) takes on values in the interval 
between s and s + ds with a probability density, V{s), that converges to a non-Gaussian 
limiting distribution, say Voo{s), that decays roughly as exp(— cs^) for large s}^ 

For the three-dimensional {d = 3) case, Bleher and Dyson^^'^"^ established that n^lQ 
grows as CVlnC when C — > oo. The behavior of s = ns{()/(y/ln( as a function of exhibits, 
just as for d = 2, essentially random fluctuations about a zero mean: For a depiction of the 



9 



behavior of the wild variation of the truncation error, -E'o(C); ior d — 3 see Figure 3 of ref 8. 
Bleher and Dyson also conjectured on the basis of a heuristic argument that the distribution 
of the variable, s — n3{Q/(^/[n(, would approach a Gaussian distribution. To illustrate this 
conjecture numerically, we present in Figure 3 the finite-C distributions, obtained for 

C = 5 up to 1000 using increments = lO"''. For ( > 100, the distribution is already very 
close to the limiting distribution which is also shown using an estimated variance — 6.85. 
To check that the distribution is indeed a Gaussian we have computed the ratio 
obtaining the values 2.52, 2.77, 2.89, and 3.00 for the ( ranges shown in Figure 3. For a 
Gaussian distribution, this ratio should be exactly 3. 

The appearance of a logarithmic factor in s (not predicted by our heuristic arguments) 
would normally be associated with an extra parameter, i.e., In^ =^ IncoC — (l^^C + Inco). 
It is not so surprising, therefore, that the data of Fig. 3 even for ( > 100 show systematic 
deviations from symmetry around s = that are consistent with a value of Cq > 1. 

Finally, we may recall that for o? > 3 the truncation errors always decrease in magnitude 
as 1/C^; and, in analogy with critical phenomena, one might speculate that the errors will 
always be distributed in Gaussian fashion. 

Ill Exact Results for One-dimensional Models 

III.l Motivation 

Here we consider exactly solvable {d — l)-dimensional fluid models in a discretized space 
and calculate the critical parameters as functions oi ( — a/ao. Recall again that a measures 
the size of particles which, in this section will always be the hard-core diameter, while Gq is 
the lattice spacing of the discretized system. We will verify exphcitly that discontinuities in 
fi{r) or its derivatives dominate the convergence of the critical parameters Tc{C) and Pc(C) 
when ^ — >• oo. 
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Now, if 2;;L) is the grand canonical partition function of a linear one- dimensional 

system of length L at temperature T = l/k^P and activity z, the appropriate one- dimensional 
"pressure" in the thermodynamic limit is given by 

(3p= lim (l/L)lnS(/5,2;L). (7) 

L— >oo 

This limit can be most efficaciously calculated via the construction of the Laplace transform 

POO 

<f{P,z;s)= / e-^^S(/3,^;L)rfL, (8) 

^0 

where the integral is absolutely convergent for all values of s with a real part exceeding 
the abscissa of convergence, So(/5, -z).^^ As s approaches Sq from above along the real axis, 
\l/(/3,2;;s) diverges to infinity as a simple pole.^*^ From eq[7|the thermodynamics are then 
fully determined by 

(3p = so{f3,z). (9) 
For a linear system with only nearest-neighbor interactions, one finds 

<il{P,z;s) = z/[l-zJis)], (10) 

where J{s) is the Laplace transform of the Boltzmann factor associated with the potential 
(p{r), namely, 

POO 

J{s)= e-''exp[-Pifi{r)]dr. (11) 
Jo 

The abscissa of convergence, Sq, is determined by the simple pole in eq^J From eqlHlthe 
equation of state may then be derived from^^ 

l/z = J{/3p), -l/pz = J\Pp), (12) 

where J'{s) denotes the first derivative of J{s) with respect to s. 

For a discretized system the integral in eq |H1 must be replaced by a sum and, similarly, 
J(s) in eqs UnHEl must be defined as 

oo 

J{s; C) = Y1 e"''"" exp[-P^{lao)]ao, (13) 
/=o 

11 



where the identification r = lao with / = 0, 1, 2, ■ ■ ■ has been adopted. 

A one-dimensional model with finite-range potentials, <f{r), such as we will consider, does 
not, of course, exhibit any phase separation or criticality. Nonetheless, questions regarding 
the critical point parameters, T^,, pc, etc., can be approached^^ by introducing an appropriate 
van der Waals limit. In this limit the equation of state may be written 

P = Po{p,T) -eap'^, (14) 

where po{p,T) is the solution of eq fT^ while the second term, eap^, is chosen so as to match 
the exact high-temperature behavior of the second virial coefficient which results when an 
infinitely weak, infinitely long-range attractive interaction is added. In the following we 
investigate various short-range interaction potentials, ip{r), and study the approach of the 
critical parameters for discretized systems to their continuum limits when ( oo. 

III. 2 Pure hard-core model 

The simplest model is the hard-core model with interaction potential 

ip{r) = oo, for r < a, 

= 0, for r > a. (15) 

In the continuum limit we obtain from eq^lthe well known result 

Ppo{p,T) = p/{l-ap), (16) 

which in eq HH yields the standard vdW equation. From this, using the normal classical 
criticality conditions 

{dp/dp)T = 0, {d^p/dp\ = 0, {d'p/dp\ ^ 0, (17) 

we find the usual results 

p*(oo) = apc(oo) = i, T*{oo) = kBTc{oo)/e = ^, (18) 
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while Pc{oo) = apc{oo)/e = ^. 

Now, consider the discretized system with parameter ( = a/aQ. Let us first address the 

case in which ( is an integer, say n (as has been standard in the hterature^"*^). It is then 

convenient to define an "effective" or lattice hard-core diameter &o = ^"^o- When ^ = n we 

have, of course, Bq = a. From eq^lwe then easily obtain 

aoPo{p,T;n) f 1 - {bo - ao)p\ , . 

— m ; , Oq = naQ = a. (19) 



From this, eqs El and El the critical parameters are found to be 

p (n) = 6oPc = -T 20 

6[n — 1) 

k T 

T*{n) = = 2apc(l - 6oPc)(l - &oPc + boPc/n). (21) 

When n — >■ cxD the continuum values in eqElare reproduced. By expanding in inverse powers 
of n, we find that the rate of convergence is described by 

^Pcjn) ^ [pc(ri) -Pc(oo)] L+ ^ (22) 

Pc(oo) pc(oo) 2n\ 4n 8^2 " 7 ' 

AT,(n) _ [Te(n) - r,(oo)] _ 1 ^ A + ^ + . . .y (23) 



Tc(oo) Tc(oo) 2n \ 8n IQ-n? 

As anticipated from the analysis of the hard-core part of the second virial coefficient presented 

in the previous section, Pc{n) and Tc{n) both converge to their limiting values as n"^, i.e., 

with exponent ip = 1. 

Apart from matters of convenience or tradition^*^ there are no reasons for restricting the 
choice of C to integers. However, one may reasonably guess that the truncation errors in the 
critical parameters will still decay as 1/C for most values of C,. But an interesting question 
arises: Are there any particular values of C, that improve the convergence rates? In other 
terms, if we put Q = n — 5 with < 5 < 1, are there 5 values such that Pc(C) and Tc(C) 
converge to their limits faster than n~^l 

To address this question, consider the second virial coefficient, B2(T), for the potential 
ip{r), in both the continuum and the lattice: see eqs El and El for d = 1. Certainly, one has 
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B^iT) = a in the continuum while B2 = {n — ^^^aQ in the discrete space when C, = n — 5. 
Note that, by definition, = a/( = a/{n — 6). It is then plausible that the choice of 6 
which matches B2 to B^ may improve the convergence of Pc(C) ^c(C)- the hard-core 
model, this value is simply 5 = |. 

Accordingly, let us repeat the previous calculations but with ( = n—^. But the hard-core 
part of the pressure, pq{p,T), given in eqUHldoes not change! Thus the same results, eqs 
1201 and EH for pc and T^, still apply. Note, however, that one now has Oq = a/{n — |) and 
= naj {n — ]^). Substituting these into eqs EOl and ED yields 

^' 'i + r + T7-^ + ---h (24) 



^TJC) 1/19 

- 'l + - + 7nr-^ + --- . (25) 



Evidently the hoped-for improvement in convergence has been realized! Indeed, not only 
does the choice 5 = \ yield faster convergence by a factor l/n^ 1/^, but the magnitudes of 
the leading amplitudes here are much smaller than those in eqs |221 and 1221 It is also worth 
noting that they are of opposite sign. Consequently we learn that convergence from above 
or from below should not be a universal feature. 

III. 3 Square-well models 

To explore further let us consider hard-core square- well models of well depth e > 0. The 
interaction potential is thus 

ip{r) = 00, for r < a, 

= — e, for a < r < c = Xa, 

= 0, for r > c, (26) 

with A > 1. Realistic values satisfy A < 2 but, for the present purposes we may also examine 
A > 2. 
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As in the previous section, let us first suppose that = n is an integer; but since A need 
not be integral it is helpful to employ the (slightly nonstandard) notation [x] = minm{m > x} 
where m is an integer, so that [x] is the smallest integer not less than x. Then the required 
Laplace transform is 

/3e-naos i ( T _ plSe\ -[\n]aos 

As-X) = ao'- , (27) 

where, once more, aQ = a/Q = a/n. The pressure po{p,T) again follows from the relations 
eq^lbut an analytical treatment is no longer tractable. Nonetheless, by solving the critical 
point conditions eq II 71 numericallv. precise values for and pc are readily found. The results 
for A = 2, 3 and 5 are presented in Figure 4. As expected from the previous analysis of 
the pure hard-core model, Pc{n) and Tc{n) converge to their limits as 1/n. 

Now the previous question again arises, namely, can one find nonintegral values of ( that 
provide enhanced rates of convergence? It is reasonable to anticipate that if such values exist, 
the corresponding second virial coefficients B2(T) [for the potential v^(r)] will also converge 
faster to the limit B^{T). Accordingly, let us again look for nonintegral values of ( = a/aQ 
for which B2(T) = B^(T). To simplify the formulae it is convenient to write c/ao = 9 and 
to suppose (as will transpire) that 6 is also nonintegral. Then we find 

BliT) = ([C]-i)ao + ([^]-[C])ao(l-e-^^), 

^ a + (c-a)(l -e"^'), when ( ^ oo. (28) 

It is possible that there are values of ( that solve the equality -B2 (^) = B^(T) and depend 
explicitly on /5 oc 1/T. However, such T-dependent choices for C are hardly acceptable when 
the actual aim is to determine the unknown value of T^. Accordingly, we will seek only 
T-independent solutions. To that end we may first let /? — 0: in this limit the second terms 
on the right hand sides of eq|2Hl vanish and the model reduces to the pure hard-core case. 
Equating the first terms in eq |2H1 yields C = [C] ~ | which is generally solved hj ( = n — ^ 
with n = 1,2,3, ■■■ . Of course, this just confirms the previous conclusion. Equating the 
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second terms (for n < oo and for n = oo) similarly yields 6 = [6] — ^ so that one must also 
have 9 = m — J with m = 1, 2, 3, ■ ■ • . 

But at this point one must observe that 

c 9 m - I 2m - 1 
a C 2n — 1 

is to be held fixed when one allows ( = n — ^ to diverge. Evidently this is possible only if 
A is the ratio of two odd integers say, in lowest terms, A = (2j + l)/{2k + 1) with j and k 
nonnegative integers and j > k (to preserve A > 1). If numerator and denominator here are 
multiplied by {21 + 1), where / is a positive integer, the odd/odd character of A is preserved. 
Comparing with eqEHlwe then quickly find that the arithmetical sequences 

C = Ckin) = {2k + l)(n - \) {n = 1, 2, 3, ■ ■ ■) (30) 

satisfy the second virial coefficient equahty provided A takes a value of the form (2j + 1)/ 
{2k + I). 

In the simplest case, /c = 0, the ratio \ = c/a can only be an odd integer. When k = 1 
further acceptable or "conforming" values are A = l|, 2|, 3|, ■ ■ ■, with ^ = 4^, 7|, 10|, 
■ ■ ■ . Likewise, k = 2 yields A = |, |, ^, ■ ■ ■ and ( = 7|, 12|, ■ ■ ■ . Figure 5 presents 
critical densities and temperatures vs for the optimal sequences C = Cfc for the models 
with A = l|, 2|, 3 and 5. The expected enhanced convergence with exponent ip = 2 is 
clearly confirmed for these special, conforming values of A. Notice, furthermore, that even 
for relatively small values of (, the magnitudes of the truncation errors in and pc are 
ten-fold smaller than for the integral ( values. Conversely, if A is not the ratio of two odd 
integers one cannot increase the convergence exponent merely by choosing some optimal 
arithmetical sequence of ( values. Nonetheless, when A is close to a conforming value the 
related (k values should yield more accurate estimate for Pc(oo) and Tc(oo). 

Finally, we must point out that despite the somewhat tortuous derivation of the conform- 
ing A values and the corresponding optimal sequences Ck{n), there really is no mystery in 
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these results! Rather one need only focus on the dominant feature of the hard-core square- 
well models, namely that the Mayer /-functions or, equivalently, the Boltzmann factors 
^-l3<fiir)^ exhibit discontinuities at r = a (as for the pure hard-core model) and at r = c: 
see Figure 6(a). The slow l/( decay of the truncation errors is simply ordained by these 
discontinuities. 

To see this clearly in the general case of approximating a single-variable integral by a 
simple sum over lattice sites spaced at intervals r/+i — = oq, suppose that the integrand has 
a discontinuity of magnitude A at r = 6. If the closest lattice sites are located at r+ = b+aQ6 
and = b — a^d' with (5 + 5') = 1 and 5,5' < 1, the leading error in approximating the 
integral through r = 6 by a sum is of order A(| — 5)ao oc Aa/C. this is easily seen graphically. 
However, if r = 6 lies midway between adjacent lattice sites, so that 5 = 5' = |, the leading 
error term, proportional to A, vanishes identically! Instead, if the integrand has bounded 
derivatives near r = b, the error is only of order Oq = jC^. 

When there are two discontinuities in the integrand (and the origin is fixed at a lattice 
site) an essentially number theoretic problem arises: namely, one must find a lattice spacing 
ao so that both points of discontinuity at r = a and r = c lie midway between adjacent 
lattice points. From this perspective the condition (^ = a/ao = n— | and the restriction 
eq EHl on A are rather obvious. Ameliorating further discontinuities would require further 
restrictions on their relative locations. 

We conclude from this discussion that if the Mayer function, /(r), is continuous (but 
perhaps with discontinuities in its derivatives), the critical parameters should converge more 
rapidly than To investigate this issue we examine first a model that can display a 

logarithmically softened repulsive core. 
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III.4 Logarithmic repulsive core models 



We define the logarithmic repulsive core potential by 



oo, 



for r < a, 




for r > c, 



ior a < r < c = Xa, 



(31) 



where the parameter w — 1 — w' controls the continuity of e~^'^^^^ ect r — a. Indeed, when 
w = the potential ip{r) diverges smoothly to +oo as r ^ a+; however, for < w < 1 the 
potential rises only to eln(l/w) < oo. Thus the Mayer /-function is continuous everywhere 
when w — while it exhibits a discontinuity at r = a whenever w ^ 0: see Figure 6(a). 
Furthermore, the derivative, df{r)/dr, is discontinuous at r = c for all values of w. One may 
then anticipate that lor w both pc{n) and Tdn) converge to their limiting values as 
while for w = the convergence should be faster. 

To confirm this, we consider for simplicity only integer values of C, and set A = 2. The 
Laplace transform of the potential in the discretized system is then 



where ao = a/n. From this we are able to calculate po(P) T) numerically, obtain the equation 
of state, and solve to obtain the critical point values illustrated in Figure 7 for different 
values of w. The critical parameters for w = clearly converge to the continuum limits 
as 1/n^, faster than the standard 1/n rate observed for the previous models with "abrupt 
hard cores" : see the insets in Figure 7. On the other hand, for w 7^ the truncation errors 
decrease, as expected, only like 1/n. 

As explained, the slow 1/n decay when w 7^ arises from the hard-core discontinuity 
of /(r) at r = a. By choosing C, = n — \ so that the hard-core part of the second virial 
coefficient for the discretized system matches the continuum counterpart, we again expect 




l=n 



(32) 
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to obtain faster convergence. Indeed, as illustrated in Figure 8, those values of ( do lead to 
a rate of convergence of Pc(C) ^c(C) to their limits. 



III. 5 Cubic model 

As a last case it is natural to consider a potential, v?(r), smoother than the previous 
examples and ask how rapidly the critical-point values approach their continuum limits: We 
expect to find a rate faster than l/n^. For this purpose, consider a cubic model defined by 
the potential 



oo, 
— eln 
0, 



r \ 
""^2 ( - - A2 



for r < a, 
for a < r < c, 
for r > c, 



(33) 



where 9, Wi, W2, and Ai, and A2 are parameters. For illustration and to ensure (p{r) and 
{d(f/dr) are continuous at the boundaries, r = a and c = Aa, we choose B = 3, A = 2, 
Ai = A2 = |, Wi = I and W2 = 2. With these choices, /(r) and {df /dr) become continuous 
everywhere: see Figure 6(b). 

The second virial coefficient for this potential with ^ = n is given by 



^2 (^) = «0 



2n-l 



n — 



(34) 



1 _ 2 ^ (e-/3^('«o) _ 1) 

l=n 

with ao = a/n. This can be calculated numerically for any temperature. It is then easily 
seen that -B2 (^) converges to B^{T) as 1/n'^: this is just as anticipated since /(r) and its 
first derivative are continuous. 

The critical point for this model can also be computed numerically via eqs IT^ and IT^ bv 
using the Laplace transform 



2n-l 



(1 - e^ne-2"<^«^ 
J{s) = ao ; h ao > e 

l=n 



■laos 



1 3 n 3 

2 ^ 2 W ~ 2 



n 2 



f3e 



(35) 



The results are presented in Figure 9. We observe that both Pc{n) and Tc(n) converge to the 
limiting values as just as indicated by the behavior of the second virial coefficient. 



19 



IV Extensions for Higher Dimensions 

So far we have discussed in explicit and analytical detail the effects of the discretization 
parameter, (, on the critical temperature and density, and Pc, of various one- dimensional 
models. A central lesson that we have learned from the exact calculations is that the behavior 
of the critical temperature and density in discretized systems closely reflects the convergence 
of the discretized second virial coefficient, -B|(T). As indicated in the Introduction and 
discussed further in Appendix A, one may extend this notion to analyze more realistic, 
higher dimensional systems, e.g., two- or three-dimensional model fluids. As we will show 
via explicit numerical calculations and simulations, the basic conclusion is confirmed. We 
aim here, in particular, to discover techniques for estimating the limiting critical points 
more precisely on the basis of less computationally expensive simulations of the discretized 
versions. 

IV. 1 Smooth potentials 

With this end in mind, let us first examine the convergence behavior of the discretized 
second virial coefficients, for typical smooth potentials ip{r) in d = 3 dimensions. It is 
convenient to normalize via 

B;{T;0 = Bi{T)/{27Tay3), (36) 

where a (= a) is the diameter of the repulsive core of the potential defined for single-well 
smooth potentials, as usual, by '^{a) — 0. Of prime interest is the Lennard- Jones (LJ) (12, 6) 
potential,^° with well-depth e, for which the original discretized simulations^ with ^ = 1, 2, 3, 
5 and 10 indicated rather satisfactory agreement with well established continuum estimates 
of Tc- However, these initial studies did not examine carefully the behavior of T^IQ as a 
function of Subsequent investigations^ employed the modified Buckingham exponential-6 
(El-6) potentiaP^'^^ and looked systematically at the dependence of Tc{C) and Pc(C) C fo^" 
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the values C = 5, 6, • • • , 10 and ( — 15. Accordingly, we will also examine B2{T; () for the 
modified E-6 potential. 

The Boltzmann factors, e~f^'^^'^\ for the two models at the estimated critical temperatures, 
k^Tc/e — 1.299 and 1.243, respectively^'^^ are presented in Figure 10. For the LJ potential 
the Mayer /-function is smooth everywhere, all its derivatives being continuous. On the 
other hand, the modified E-6 potential exhibits a discontinuity at rmax — 0.2285(7 owing to 
the hard-core cutoff. However, the discontinuity is invisible in this figure and irrelevant in 
the computations since its magnitude is less than 10~^^^^. 

Figure 11 presents the reduced, discretized second virial coefficients for the LJ and mod- 
ified E-6 potentials evaluated for integral values of C < 25 at the best estimates of the 
corresponding critical points. Since, of course, {T; () is an analytic function of T for 
all C < oo, the precise value of T should have a rather small infiuence on the C variation. 
As anticipated from the smoothness of the LJ and modified E-6 potentials (see Figure 10), 
both second virial coefficients converge rapidly to their continuum limits. (The latter are 
explicitly marked on the figures; needless to say, the continuum second virial coefficients for 
such radially symmetric potentials can be computed to high accuracy using standard inte- 
gration routines.) This figure demonstrates that over a computationally relevant range of 
discretization parameters (i.e., ( < 20) the observed rates of convergence are well described 
by a dependence. Indeed, as mentioned in the Introduction, an effective convergence 
exponent ~ 6 has been observed numerically in the variation of Tc{C) for the modified 
E-6 potential.^ (For this model, to within the precision of the simulations, no significant 
variation of Pc(C) could be resolved.^) 

Although the agreement between the observed convergence behavior of Tc(C) and the 
second virial coefficient, i?2, in a more-or-less realistic three-dimensional model is gratify- 
ing, the value ^ ~ 6 should not be taken too seriously! Thus, it would be interesting to 
see how far the same effective rate of convergence {t/j ~ 6) is maintained for both models 
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when the evaluation of the discretized second virial coefficients for the LJ and modified E-6 
potentials is pushed further (and, hence, beyond any normally relevant level of precision). 
However, significant further efforts are needed to resolve the issue numerically and we have 
not undertaken the task. Of course, for extremely large values of ( we expect a ^/^nX/C'^ 
asymptotic convergence law for the modified E-6 potential to appear in view of the sharp, 
although exceedingly small discontinuity in /(r). 

IV. 2 Potentials with hard cores 

Because of their conceptual simplicity and theoretical interest, it is important to address 
models where ip{r) consists of a hard-core potential (fo{r) plus a bounded and smooth at- 
tractive potential <fi{r). Such systems include both the RPM and the hard-core dumbbell 
models studied in ref 5. Certainly the convergence of the second virial coefficients will be 
dominated by the behavior of the hard-core part Bq (which, of course, is T- independent) . 
Thus, following the arguments for d = 1, we may expect that the critical temperatures and 
densities will approach their continuum limits as ^ — > oo in a manner reflecting the variation 
of Bi 

As seen in section II the hard-core truncation error Eo(T) = {Bq/B^) — 1 decays as 
for all (i > 4 while for d = 3 the error varies as (In^)^/^/^^. However, as demonstrated 
in Figure 2(b) (for d — 2) and as shown in Figure 3 of ref 8 (for d — ?>), the behavior 
of the truncation error for d > 2 is intrinsically erratic, chaotic or 'noisy'. Indeed, in 
d — ?) dimensions it may be viewed as a random variable with, as discussed, a Gaussian 
distribution: recall Figure 3. Furthermore, it is found in the simulations that the estimates 
of Tc(C) and Pc(C) c.g.,^ the RPM, also exhibits noisy behavior that, at least superficially, 
is reminiscent of that which necessarily arises in Bq. 

To proceed we might initially ask, as in section HI, if one might not be able to select 
special sequences of values for C, (not necessarily integers) so that T^^Q and Pc(C) would 



22 



converge faster than otherwise. To that end one should first ask (in d > 2 dimensions) for 
sequences of ( values for which = so that Eq{Q vanishes identically. Such matching 
sequences would seem good candidates for ensuring the optimally rapid convergence of Tc{Q 
and Pc(C) 3S verified in section III ior d — 1: see Figures 5 and 8. 

Inspection of Figure 2(b) and Figure 3 in ref 8 leaves little doubt that such infinite 
matching sequences exist ior d — 2 and 3. But one is easily convinced that integral values 
of C cannot appear; further, one may seriously doubt that even arithmetical sequences like 
kn — S with fixed S and integer k could satisfy the matching condition. However, there is no 
obstacle to the computation of matching values to any required precision up to relevantly 
large values of C- Explicitly, in Table I we present, for both d = 2 and 3, selected matching 
values of ( in the range (5, 25). In planning simulations of models with hard cores it seems 
likely to be advantageous to choose such matching values. However we have not undertaken 
simulations to verify this surmise. 

It maybe worth mentioning that since periodic boundary conditions^° for a range of box 
sizes must be employed (in order that essential finite-size extrapolation techniques can be 
applied systematically^^) a purely cosmetic blemish will arise if matching values of ( are 
adopted: specifically, since L*^ must be an integer (where L* = L/a is the reduced linear 
dimension of the simulation box),^° nonintegral and, indeed, irrational values of L* will be 
entailed. However, this feature is of no consequence in any of the computations entailed 
in implementing the fine-discretization technique,^"^ neither setting-up, data gathering nor 
analysis. 

However, whether or not matching C- values are employed for hard-core systems (and 
whether or not they prove efficacious to some degree in practice) it is reasonable to ask 
how the difficulties of extrapolation 'through inherent residual truncation noise' might be 
ameliorated. Explicitly, faced with evidently erratic estimates for T'c(C) Pc(C) such as 
found for the RPM in ref 8 (or for hard-core dumbbell models in ref 5), how might reliable 
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estimates for the continuum parameters be obtained? 

Following ref 8 we describe two approaches. First, noting that the general behavior of 
Tc(C) and Pc(C) resembles that of Bq and accepting the idea developed in section III that 
this should be more than coincidental, let us attempt matching post facto. Thus we may 
introduce a modified or smoothed discretization level, CKO defined so that, for d = 3, 

Eo{0=c'^/C\ (37) 

where c'^ is an assigned constant. (Here we neglect, for the time being, the factor (In^)^/^ 
known to be present in the exact asymptotic behavior.) By construction, the truncation 
error -E'o(C) iiow decays smoothly as We may thus hope that T'c(C) and Pc(C) will 

also behave smoothly when plotted vs 1/[C^(C)]^. In order for this method to be convincing, 
however, the resulting extrapolated estimates should not depend significantly on the choice 
of c^. To demonstrate that this is indeed so for the RPM, we present in Figure 12 plots of 
Tc(C) and Pc(C) as functions of C^(C) computed with the choices = 15i?o(5), 25-E'o(5), and 
35^0 (5). [We remark that in ref 8 only = 25Eo(5) was used.] Evidently the behavior is 
fairly smooth and, in fact, much less erratic than plots vs. see Figure 4(i) in ref 8. From 

these plots we conclude T*{oo) = kBT^Da/ql = 0.04933(5) and p*(oo) = p^a^ = 0.0750(10) 
for the continuum RPM.* 

Of course, the presence of the logarithmic factor in the true asymptotic decay of -E'o(C) 
raises a question regarding the validity of this smoothing procedure. However, one may 
anticipate that the logarithmic factor plays an insignificant role for the fairly small values 
of C = 5-20 employed in these simulations of the RPM. Indeed, one may readily introduce a 
factor (In^)^/^ in the definition of C^(C) eqEZl but on doing so we find that to within the 
simulation uncertainties the approach provides the same estimates for Tc(oo) and pc(oo). It 
goes without saying, however, that the effectiveness of this method for the RPM is not much 
more than a promise of success in other cases; but reasonable optimism seems in order! 
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IV. 3 Rescaling to accelerate convergence 

Another approach to improving convergence that we have explored to a hmited de- 
gree is based on the idea that what is needed in obtaining optimal estimates from a fine- 
discretization calculation is some appropriate rescaling of density and temperature depending 
on (. Our analysis again suggests that this might profitably be guided by the behavior of 
the full discretized second virial coefficient, 



rather than just the repulsive part Bq(T). Indeed, the contribution of the attractive tail fi{r) 
must surely play a role and a hard core might be softened as in the LJ and E-6 potentials. It 
is then natural to examine the behavior of the rescaled discretized critical densities defined 
via 



in the hope that they will vary more smoothly vs l/C^ for appropriate tp. 

To go a little further, one might, in as far as (pi{r) varies relatively smoothly by definition, 
approximate the full second virial coefficients in eqEHl by replacing the contribution B'i{T) 
by a ^-independent value, say bi. Indeed, when one considers the RPM, some modification 
of eq EHl is essential since, by virtue of the long-range Coulomb interactions, the first low 
density correction to ideal gas behavior is universal and proportional to p'^/^ rather than 
to p^. And although a term also exists, it is dominated by a In p term.^^ Naively, one 
might be tempted to say that the second virial coefficient is divergent; but in reality although 
a hard-core term of the form Bq{T) plays a part, there is no direct meaning that can be 
ascribed to Bi(T). However, if B^iT) in eqEHlis replaced by a parameter b1, we may regard 
the density rescaling eq EHl merely as a semiphenomenological device. A similar expression, 
but with a parameter b^, may then be used to rescale Tc{() to obtain a T^{C). Then the aim 
should be to find suitable choices of 6^ and 6f that provide enhanced convergence. 



B^,{T) = B^,{T) + Bl{T), 



(38) 



pt(C)=Pc(C)s^(T)/5r(T) 



(39) 
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This procedure has been successfully demonstrated for the RPM in ref 8. There it is seen 
(in Figure 4) that the choices 6^ = 0.45^ and bj = 1.7B^ yield plots of pJ(C) and Tj"(C) 
that display very little ( dependence in the range ( — 5-20 (whether plotted vs. as 
seems most appropriate, or employing some other ip values). Indeed, extrapolation of pJ(C) 
and T^{C) is then limited almost entirely by the Monte Carlo simulation uncertainties at the 
given C values: and the same estimates are obtained for the RPM as quoted above. 

Finally, as the previous remark brings out, the practical application of the extrapolation 
methods just outlined to simulations of model fluids relies on the use of methods capable 
of providing precise and reliable estimates of Tc and pc at a number of fixed discretization 
levels. On the other hand, although the approach is as yet untested for = 2 and 3, the use 
of the favored matching discretization levels listed in Table I may provide the best strategy 
if an investment in only one or two simulations at chosen ^ values can be justified. 

V Summary 

In conclusion we have studied the effects of the lattice discretization level, C, on the crit- 
ical parameters, in particular, the critical temperature, Tc, and density, Pc, of model fluids. 
The lattice discretization technique can elevate the speed of simulations for complex fluids by 
factors of 10-100 relative to the corresponding continuum models. It has been demonstrated 
that the convergence of Tc{C) and Pc(C) as ( increases may be understood rather straight- 
forwardly by studying the discretized second virial coefficient, -B^(T). Appendix A discusses 
this conclusion from a general theoretical perspective. In particular, when the interaction 
potential contains a hard-core repulsive part, the partial second virial coefficient, B^, for the 
hard-core contribution (which is essentially given by the number of lattice points inside a 
d-dimensional sphere of radius () dominates the convergence of i?|(T). Rigorous aspects of 
this issue, namely how the number of lattice points converges to the volume of the sphere, 
were presented in section II. The relative truncation error in -62(7), namely £'o(C) as defined 
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in eq was shown to decay as q/C*-'^^^''/^ for d < 3 with Cd ~ 0(1) for (i = 1 and 2 but 
C3 ~ A/ln~C, while for all > 4 one has Eq{() ~ These results can be understood 

intuitively on the basis of heuristic arguments presented in Appendix B. 

To illustrate the general arguments in further detail, we have studied various one-dimensio- 
nal model fluids of van der Waals character in section III obtaining exact results for the 
discretized critical points. Indeed, it was confirmed rather generally that the ^-dependence 
of Tc{C) and Pc(C) reflects closely the behavior of the second virial coefficient. Further- 
more, the smoothness of the potential, f{r), or more precisely of the Mayer function, 
/(r) = e~'^^^^^''^'^ — 1, governs the improved rate of convergence of the critical parame- 
ters when /(r) is continuous. Thus d = 1 systems with a hard-core potential and jump in 
/(r), generally exhibit only a slow, convergence of Tc{() and Pc(C) to their —>■ 00 limits. 
However, if ( is chosen properly — in general not as an integer — to match the second virial 
coefficient, B2, to its continuum limit, B^, then the convergence becomes faster by a factor 
Furthermore, the relative deviations of the critical parameters from the continuum Tc 
and pc values become much smaller: see sections III.B-D. When /(r) and its derivatives are 
continuous, the convergence rate is improved by further factors of 1/^: see sections III.D 
and E. 

The exact results for the one- dimensional models provide a useful guide to what may be 
expected in higher dimensions. However, understanding the behavior of more realistic two- 
and three-dimensional models requires further considerations that are presented in section 
IV and confirmed numerically via various simulation studies. From the results of section 
II, one learns that the second virial coefficients for models with hard-core potentials vary 
in an intrinsically erratic fashion as ( increases. Thus even though the overall magnitude 
of the truncation error, Eq((), decays as (except for a logarithmic factor when 

d = 3), its wildly erratic behavior severely hampers extrapolation of the critical parameters 
to their continuum limits for model fluids with hard-core potentials. To improve matters, a 
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smoothed discretization parameter, C^(C)) was introduced that serves to smooth the erratic 
critical data in a controlled way and thereby ease extrapolation to the continuum limit. A 
flexible procedure for resettling p and T on the basis of Bq was also demonstrated successfully 
using simulations for the restricted primitive model electrolyte. 

On the other hand, it was shown that it is feasible to select special values for ( that 
make Eq{() vanish: see Table 1. By adopting such optimal values one may hope to obtain 
critical point estimates lying sufficiently close to the continuum values that the need for 
extrapolation to C = C)0 will be obviated in many practical cases. 
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A Relating Virial Coefficients to Critical Parameters 

If eqEl the rather general but approximate equation of state of van der Waals character, 
is accepted as a reliable guide, a direct proportionality between the rate of convergence of the 
discretized second virial coefficients, as defined in eqlH and the critical parameters, Tc(C), 
Pc(C); etc., must be expected as the discretization level, (, approaches infinity. 

However, this conclusion has a much broader basis. To see this, consider the Mayer 
fugacity expansion 



oo 




(Al) 



1=2 
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or the virial expansion 

oo 

p/k^T = p + Y,Bi{T)p\ (A2) 

1=2 

and recall^^ that for a system with pair interaction potential ^{r) the coefficients hi{T) and 
Bi{T) can be defined exphcitly as multiple integrals over products of Boltzmann factors 
exp[— (/9(r.y)/A;B^]- Of course, these results hold in all dimensions d and can be extended 
to multicomponent systems, many-body interactions, etc. In a discretized system all the 
cluster integrals must be replaced by sums as in eqlH 

The series in eqs lAll and IA2I and their discretized versions converge at small z and p for 
models of practical interest and in principle, by analytic continuation, then yield Tc and pc- 
For example, for systems displaying nonclassical critical behavior, which are of prime interest 
here, it suffices to examine the inverse compressibility, 1/Kt{T, p), which vanishes uniquely 
at T = Tc and p = pc- (Spinodal divergences of Kt do not arise in realistic models owing to 
the essential singularities encountered at points of condensation^^ where the compressibility 
remains finite.) 

As C increases the discretized second virial coefficient B2iT) will approach its limit, say, 
with an exponent tpB- Unless accidental cancelations occur, one must, via eqs lAll and IA21 
expect that the pressure and all thermodynamic quantities derived from p{T, z) or p(T; p) in 
the usual ways by differentiation, etc., can converge when C, ^ oo no more rapidly than do 
the additive terms h2{T)z and B2(T)p [where, of course,^'' -B2 = —^2 ]• Conversely, unless 
the Boltzmann factors for any many-body potentials, which enter the expansions at higher 
order in z and p, vary more sharply, e.g., with derivative discontinuities of lower order than in 
exp[—(f{r)/kBT], the higher-order multiple discretized cluster integrals should not introduce 
any slower convergence factors. Indeed, further integration will generally have a smoothing 
effect leading to enhanced rates of convergence: see the dependence on dimensionality d 
discussed in Appendix B. We thus conclude that, in general, the convergence exponents ipT 
and ipp match ipB (or in very special cases — see the discussion related to Table I — may 
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exceed ips)- 

These arguments might seem to side-step the issue of nonclassical criticahty raised in the 
Introduction. The crucial point, however, is that the second virial coefficient (discretized or 
not) closely reflects and, essentially, embodies the interaction potential (p{r). But, as found 
explicitly even in highly nontrivial exactly solved models with nonclassical critical behavior 
— in either one dimension^^ or in higher dimensions'^ — the critical temperature Tc varies 
smoothly and analytically with the potentials. Likewise for pc when this is not trivially 
fixed by some symmetry, etc. A well known example, the plane triangular Ising modeP® 
with three, positive, negative or vanishing exchange couplings, Ji, J2 and J3, illustrates this 
well and also indicates how the assertion may fail in special cases! Specifically, for certain 
particular values of the potentials the nature of the critical point itself may abruptly change. 
This happens in the triangular Ising model with Ji, J2, J3 < where an anomalous critical 
point arises atT = when Ji = J2 = Js- In general, such points are multicritical points?^ 

If a discretization process is performed on a system in the close neighborhood of any 
multicritical point, one must anticipate that the C-dependence of the second virial coefficient, 
S|(T), may be an unreliable guide to the behavior of the critical parameters. Indeed, near a 
bicritical point, for example, there are critical loci, say Tc{g), which vary as \g — gy^l'^, where 
gi, is the bicritical value of a potential parameter, g, while 4> is a nonintegral (and nontrivial) 
crossover exponent. Of course, in the situations addressed in this article, the implicit 
assumption — that we may highlight here — is that the critical point under consideration 
is, as typical in applications, one that resides in a manifold of critical points with uniform 
(and universal) critical behavior. 

The insight into the nature of multicriticality reveals the fallacy in what might other- 
wise appear to be a cogent counter-example to our basic conclusion. Suppose, recognizing 
explicitly the functional dependence on the interaction potential as embodied in the second 
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virial coefficient, one regards the pressure as given by 



p = V{p,T-B',{T)). 



(A3) 



Then for nonclassical criticahty, the dependence on the ffist two arguments at (pc^c) must 
be nonanalytic; but one might suppose that the third argument also enters in some singular 
fashion. Then, from the vanishing of 1/ Kt{T, p) or by eliminating p between the standard 
conditions in eq 17, one might expect to derive a relation for Tc{C) of the form 



Now if }C{x; y) is analytic in both arguments at and through x = Xc = and y = yc = 
B^[T^), our previous conclusion that the behavior of [Tc(C)-T~] reflects [5^(Tc) -^^^(Tc)] 
is conffimed. Conversely, if /C(x; y) is singular in y at yc, perhaps varying as {y — yc)^ with 
X 7^ 1, the conclusion would generally be false: the convergence exponent tpT for would 
rather be expected to depend on some combination of ipB and x- However, the defect of this 
argument is that it presupposes, in essence, that the critical point in question at (T^,p'^) 
or {xc,yc) is in reality a special critical point sensitive to the specific potential that yields 
B2{T^) = yc- As explained, this violates our basic presupposition that, in truth, must be 
required of any general scheme of analysis such as concerns us. 

While we believe our arguments are convincing as regards the leading asymptotic depen- 
dence of Tc{() and Pc(C) matching that of B2(T), one may still ask whether the higher order 
corrections when ( ^ oo might not depend in some way on the bulk critical exponents. It 
is not implausible that the answer is sensitive to the details as to how precisely one inter- 
prets the argument of -B|: i.e., whether as supposed in eq IA4I or in some other reasonable 
way. A renormalization group approach and general scaling ansatz (such as employed in the 
articles by Kim and Fisher cited in ref 11) should be able to cast light on this issue. But, in 
the absence of some concrete and practical application, the point seems too academic to be 
worth pursuing at this juncture. 



/C[Te(C);i?K^c(C))]=0. 



(A4) 
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B Approximate Estimation of the Hard-core Trunca- 
tion Error 



Consider a c?-dimensional sphere of radius R = (a^ centered at the origin of a discretized 
space with a lattice spacing oq. Each lattice site may be regarded as the center of an 
elementary cell of volume Oq. The surface area of the sphere is = CdR'^~^ where Cdd = Vd 
is the volume of a unit sphere. This surface then intersects the elementary lattice cells and 
divides them into two parts. The elements of the surface which are cut by individual lattice 
cells have a mean area, Aq ~ g^ao"^ for C = i?/ao ^ 1 where qd is a real number which we 
expect to remain of order unity for all (. The mean number of such surface elements is thus 
given by 

Sr ~ An/Ao = {Cd/qdX'-'. (Bl) 

For d > 1 each surface element cuts a lattice cell at a different mean tangential angle 
and, depending on the location of the individual cell, its center will lie either just inside or 
just outside the sphere. When, overall, the surface elements are uncorrelated around the 
sphere, one may expect the probability of the center of a surface cell lying in each region to 
approach ^ as ( oo. Consequently, each more or less independent surface cell contributes 
a deviation or "error" of +| or — | to the total number of enclosed lattice points, Nd{C)- 
The total mean square error should then be proportional to the number of the lattice cells 
cut by the surface of the sphere and so given by 

(AiV,2) = {nl) ^ Sn6l (B2) 

where ANd{() = Nd{C) — Vd{() = nd{C)- The root-mean-square deviation 6d should be of 
rough magnitude | although correlations between nearby cells should increase the actual 
value. On using eq IBH we thus conclude that the root mean square error should vary as 

= \[^) ^ SdVC^^C^'-'^^'- (B3) 
32 



Finally, this implies that the truncation error, -E'o(C) = ^rf(C)/^rf(C) ~ should decay as 

As remarked, this result is (rather trivially) exact for d = 1 while for d = 2 it is confirmed 
by Bleher et al}^ For d = 3, however, the exact results^^ imply that an additional logarithmic 
factor, ln(^, should appear in eg IB21 Furthermore, the argument fails for d > 3 since n(i{C,) 
grows like C'^^^ implying^^'-*^^ -£^0(0 ~ 

To understand why there should be a borderline or crossover dimension we may, first, note 
that if, in d dimensions we had considered oriented hypercubic hard cores of sides 2R (rather 
than spheres) a very different result would pertain. Indeed, a d-hypercube can be decomposed 
into 2d rectangular hypercones of opening angle 6 = tc/A about each lattice axis. The face 
(or "base") of each pyramidal hypercone would be of area = {2RY~^ = (2ao)°'~^C°'~^ ^ind 
on increasing R by one lattice spacing oq, the number of enclosed lattice sites would increase 
discontinuously in a completely correlated manner by, in total, 

AiV,(C) ^ 2rf(A,/a^i) = 2^rf C'^^ (B4) 

On dividing by the hypercube volume = 2'^R'^, the rms truncation error is seen to be 
Eq{() (X d/( and so decays like l/( for all d. 

This transparent result suggests that for sufficiently large d the growing surface of a 
hypersphere may, from the perspectives of the 2d outgoing lattice axes, resemble a somewhat 
flat base of a cone of opening angle, say, 9 < 7r/4. Let us, to be specific, focus on the positive 
X axis and consider the lattice layer at or closest to the (normal) distance X = R cos 6 from 
the sphere center and of transverse dimension R± = R sin 6. Then if > 3 and X increases by 
a lattice spacing ao, a fully correlated set of sites contained in a "ring" or torus of thickness 
AX = ao, width AY ~ aocot^ (provided 9"^ > 1/C) and "rectangular circumference" of 
magnitude {2R±Y~'^ is added to Nd{(). Thus fiuctuations in Nd{C) of magnitude at least 

ANj>{C) oc alR^-\sm fff-^ cos ^/a^ oc C^-^ (B5) 
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must be anticipated. But, clearly, this implies that Ed{C) cannot decrease faster than 
for d > 3, just as the rigorous results establish! Incidently, formally maximizing AA*"^ on 6 
yields 6 = 7r/4 for d > 5, 6 ^ sin-^(3-V2) ~ 35.3° for d = 4, and 6 ^ ^^nin > for 
d = 3. 
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Table 1: Selected matching values of C = o/oo in the range (5, 25) (recorded to ten decimal 
places) for d = 2 and 3 dimensions. Using these vahies in simulations of fluids with hard- 
cores of diameter a, should enhance the degree of approximation of the continuum model by 
the discretized version. 
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Figure Captions 



Figure 1: Schematic illustration of the fine- lattice discretization process for C = 1, 2, 3, and 
5 for a. {d = 2) -dimensional system where a is the hard-core diameter while ao is the lattice 
spacing. The dashed circle of radius a = (gq encloses the sites, shown shaded, excluded 
by hard-core repulsion from occupation by neighboring particles. In d — 3 dimensions these 
numbers are are N^iQ = 1, 27, 93, 251 and 485 for C = 1 to 5. 

Figure 2: Relative hard-core truncation errors, -E'o(C), scaled by for (a) a one- 

dimensional discrete system and (b) a o? = 2 simple square lattice. 

Figure 3: The distribution, V{s), of the fluctuating variable, s = [Nd{C) — Vd{C)]/Cy/ln (, for 
d — 3 calculated numerically from ^ = 5 up to 1000 using an increment = 10~^ and a 
bin size As = 0.1. Even for ( > 10 the distribution is very close to the hmiting Gaussian 
distribution that is labeled ^ = oo. 

Figure 4: Plots of (a) Apc{n) / pc{oo) and (b) ATc(n)/Tc(oo) vs 1/n for one- dimensional 
square- well models of range c — Xa with A = l|, 2, 3 and 5 where the discretization 
parameter C takes integer values. 

Figure 5: Plots of (a) Apc(C)/Pc(oo) and (b) ATc(C)/Tc(oo) vs 1/C^ for one-dimensional 
square- well models with A = l|, 2|, 3 and 5 where the discretization parameter runs through 
the nonintegral values ( — n— ^iorX — 3 and 5 but only ( — 3{n — ^) for A = l| and 2|. 

Figure 6: The Boltzmann factor, e~f^'^^'^\ vs distance for (a) one-dimensional logarithmic 
hard-core models (sohd lines) with A = 2 and w — 0, |, | and 2 at a temperature set by 
Pe — 1; the dotted lines are for the hard-core square- well models with A = l|, 2 and 3. (b) 
Boltzmann factor for a cubic model at /3e = 1 with the parameters given in the text. 
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Figure 7: Plots of (a) Apc{n)/ pc{oo) and (b) ATc{n) /Tc{oo) vs 1/n for one- dimensional 
logarithmic hard-core models with A = 2 and w = 0, 0.75, and 2 when the discretization 
parameter ( takes integer values. The insets show the behavior of the critical temperature 
and density when w = as functions of 1/rP: see Figure 6. 

Figure 8: Plots of (a) Apc(C)/Pc(oo) and (b) ATc(C)/Tc(oo) vs 1/C for the d = 1 logarithmic 
repulsive-core models with A = 2 and w = 0.5, 0.75 and 2 for the choice ( = n — ^ with 
n = 5, 6, • • • . 

Figure 9: Plots of (a) Apc{n) / pc{oo) and (b) ATc{n) /Tc{oo) vs for a one-dimensional 
"cubic" model (see text) where ( takes successive integer values demonstrating rapid con- 
vergence for smooth potentials in ci = 1 dimensions. 

Figure 10: The Boltzmann factor, e~^'^^^\ vs distance for the Lennard- Jones (12,6) poten- 
tial (solid curve) and for the modified Buckingham exponential-6 potential (dashed curve) 
evaluated at the estimated critical temperatures, k^Tcjt ~ 1.299 and 1.243, respectively. '^'^^ 

Figure 11: The reduced second virial coefficients, -82(^0; C) vs 1/C^ for (a) the (12, 6) Lennard- 
Jones (LJ) and (b) modified Buckingham exponential-6 (E-6) potentials in c? = 3 dimensions. 
The values have been calculated at the estimated continuum critical temperatures (see Figure 
10) for 5 < C ^ 25. In simulations using the E-6 potential estimates of Tc{C) have been 
observed^ to approach the limit as 1/^^^^. The insets show the same data at a larger scale 
vs 1/C^ (see also the text). 

Figure 12: Plots of (a) T*{C) and p*(C) for the RPM electrolyte vs l/C^ where C^(C) is the 
smoothed discretization parameter calculated from the hard-core second virial coefficient: 
see eq 1371 and the following text. The symbols from the left correspond to the choices 
= 35i?o(5), 25i?o(5), and 15-E'o(5). The dashed lines are to guide extrapolation to C = cxd; 
the data come from ref 8. 
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Fig. 5 
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Fig. 6 
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Fig. 8 
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Fig. 9 
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Fig. 10 
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